Time-dependent Gutzwiller theory of pairing fluctuations in the Hubbard model 
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We present a method to compute pairing fluctuations on top of the Gutzwiller approximation 
(GA). Our investigations are based on a charge-rotational invariant GA energy functional which 
is expanded up to second order in the pair fluctuations. Equations of motion for the fluctuations 
lead to a renormalized ladder type approximation. Both spectral functions and corrections to static 
quantities, like the ground-state energy, are computed. The quality of the method is examined for 
the single-band Hubbard model where we compare the dynamical pairing correlations for s- and 
d-wave symmetries with exact diagonalizations and find a significant improvement with respect to 
analogous calculations done within the standard Hartree-Fock ladder approximation. The technique 
has potential applications in the theory of Auger spectroscopy, superconductivity, and cold atom 
physics. 
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•4— > 

CO 



I 

C 

O 

o 



in 
o 

O 

oo 
o 



I. INTRODUCTION 

The present interest in the physics of strongly corre- 
lated fermion systems is accompanied by a 'revival' of 
the Gutzwiller wave function (GWF)i in order to study 
Hubbard-type models. In these Hamiltonians, doubly 
occupied sites contribute with an energy U to the total 
energy, though their weight is partially reduced by the 
Gutzwiller projector. Originally, the projector was ap- 
plied to a Slater determinant describing a spatially uni- 
form Fermi liquid. Along the years, states with long- 
range order have been considered, including projected 
BCS wave functions JLiLSik used in the context of high-T c . 

Analytic evaluations of expectation values using the 
GWF are only possible in one 6 - and infinite dimensions AS 
In the latter limit, one recovers the so-called Gutzwiller 
approximation (GA) first introduced by Gutzwiller him- 
self in Ref. ,9f. The GA, later on rederived as a saddle- 
point within a slave-boson formulation of the Hubbard 
model^ yields an energy functional for quasiparticles 
with renormalized hopping amplitude. Therefore, it of- 
fers a simple and intuitive picture for the interplay be- 
tween local correlation and electron kinetics and it is used 
in a variety of fields including the description of inho- 
mogeneous states in cupratesjii band structure calcula- 
tions^ and the theory of He 3 * 1 ^ 

In the past few years, two of us have developed a 
scheme which allows the computation of Gaussian fluc- 
tuations on top of the GA i 14 : 15 ' 16 Within this so-called 
time-dependent GA (TDGA), it is possible to evalu- 
ate dynamical correlation functions in the charge and 
spin channel, which over a wide parameter regime are 
in good agreement with exact diagonalization results 
and constitute a significant improvement over the tradi- 
tional Hartree-Fock plus Random-phase approximation 
(HF+RPA). The TDGA has been used in order to com- 
pute dynamical properties of inhomogeneous states in 



cuprates. Results obtained in this way for the optical 
conductivity^ and magnetic susceptibilit y 18 : 19 of stripe 
ordered phases have turned out to be in excellent agree- 
ment with experimental data. 

In the present paper, we generalize the TDGA towards 
the inclusion of pairing fluctuations. The significance 
of such correlations is probably most prominent in the 
context of superconductivity where the appearance of a 
singularity in the pair susceptibility signals the onset of 
Cooper-pair condensation. Here, we aim to study the 
spectrum of pairing correlations in the normal state of 
the Hubbard model, an issue which has also been ad- 
dressed, among others, within exact diagonalizatio n 20 : 21 
and various Monte Carlo methodsiSSiS&SiSiiSSiSLSS In 
the particle-hole channel, the coupling to an electric or 
magnetic field yields the optical conductivity and mag- 
netic susceptibility, respectively. In the particle-particle 
channel, direct measurements of pairing correlations are 
not so common since they do not couple to a classi- 
cal field. However, principles for the measurement of 
the pair susceptibility in metals have been discussed in 
Refs.iHiillil 

The basic idea is to probe the fluctuat- 
ing pair field of a metal in the normal state which is cou- 
pled to a superconductor via the tunnel current-voltage 
characteristics. In addition, the pairing correlation func- 
tion for local pairs is the main ingredient in the theory 
of Auger spectroscop y 33 : 34 ' 35 : 36 and also has relevance in 
the field of ultra-cold atom physics^ 

In materials which can be described within Hubbard- 
type models, a strong local repulsion induces so-called an- 
tibound states (also known as two-particle resonances), 
above the two-particle continuum. The resonances ap- 
pear as atomic-like features in the Auger spectrum. In 
Ref. [H we have shown that, despite its computational 
simplicity, the TDGA yields an excellent description of 
the two-particle response. In particular, it describes well 
the energy gap between band and antibound states and 
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the relative spectral weight even far from the dilute limit 
in contrast to the 'bare ladder approximation' (BLA), 
which is restricted to the low-density regime . 39 ' 40 Thus, 
the TDGA allows one to extend ladder-type theories 
much beyond their supposed limit of validity. The re- 
sult is an effective ladder theory where quasiparticles 
get heavier, as usual due to correlations, and at the 
same time the effective interactions between quasipar- 
ticles become strongly renormalized. These vertex and 
self-energy corrections are consistent with each other and 
do not suffer from the pitfalls frequently found in dia- 
grammatic computations, where an improvement at the 
level of the self-energy alone leads to a degradation of 
the overall performance of the theory due to the lack of 
appropriate vertex corrections ! 35 ' 36 

In this paper, besides a thorough derivation of the 
'pairing TDGA', we extend the approach to intersite cor- 
relations with extended s-wave and d-wave symmetry. 
We will show that the TDGA yields an effective inter- 
action between quasiparticles which is renormalized with 
respect to the bare U due to correlations but it does 
not include enough fluctuation effects to induce a super- 
conducting instability. This is because we start from a 
paramagnetic state treated in the GA which does not 
have enough variational freedom to describe the scale J 
of magnetic fluctuations. In addition, within the RPA 
treatment for a paramagnet the particle-particle and 
particle-hole channel are decoupled and do not influence 
each other. For spin-density wave ground states, the 
TDGA induces attractive interactions between nearest- 
neighbor pairs which are not present in the HF approach 
based on BLA. However, these are not enough to pro- 
duce a superconducting instability. Superconductivity 
due to an electronic mechanism requires the feedback of 
particle-hole fluctuations on the particle-particle chan- 
nel and this goes beyond our approach. On the other 
hand, at the level of spectroscopic quantities this is a mi- 
nor effect and our approach turns out to be in excellent 
agreement with exact diagonalizations. Here and below 
we use the terms "ladder-type fluctuations" , "particle- 
particle RPA fluctuations" and "pairing fluctuations" as 
synonymous. 

This paper is organized as follows: The formalism is 
presented in detail in Sec. HT1 where as a first step we de- 
rive the charge-rotational invariant GA functional from 
the Hubbard model in Sec. Ill Al Based on this functional, 
we show in Sec. Ill Bl how ladder- type fluctuations can be 
incorporated into the approach and how dynamical pair 
correlations can be computed. Results are presented in 
Sec. IIIII where we first exemplify the method by means 
of a two-site model. Then the dynamical pairing corre- 
lations for different symmetries are computed on small 
clusters and compared with exact diagonalization results 
and HF+BLA computations. Finally, we conclude our 
investigations in Sec. IIVI 



II. FORMALISM 
A. Charge-rotationally invariant GA 

The starting point is the one-band Hubbard model: 

where Ci a (cJ CT ) destroys (creates) an electron with spin a 
at site i, and hi a — c\ a Ci a . U is the on-site Hubbard re- 
pulsion, tij denotes the hopping parameter between sites 
i and j and /i is the chemical potential. 

The Gutzwiller variational wave function can be writ- 
ten as 

|$ G )^IP^>' ( 2 ) 

i 

where Pi partially projects out a doubly occupied state 
at site i from the uncorrelated wave function \<j>). In the 
traditional Gutzwiller approach^ the latter is a Slater 
determinant and the associated density matrix only con- 
tains the normal part: 

Here, we will consider a more general formulation in 
which is a Bogoliubov vacuum and define the anoma- 
lous part of the density matrix 

K u = I c j l^t !</>)■ 

In the general case the normal part can describe charge- 
density wave and spin-density wave broken symmetries. 
We will allow the ground state to have these broken sym- 
metries but we will assume it is normal so our saddle 
point anomalous density matrix will vanish. We denote 
quantities at the saddle point by a 0, thus = 0. The 
anomalous part is important in order to obtain the fluctu- 
ations. Indeed, in the following, we consider an external 
time-dependent perturbation which induces pairing fluc- 
tuations on \4>), then the instantaneous \<j)) will take a 
BCS-like form. 

The charge-rotationally invariant Gutzwiller func- 
tional for general charge and spin textures is derived in 
the Appendix [X] exploiting the well known equivalence 
between the slave-boson method and the Gutzwiller ap- 
proach and following previous works i 41 ' 42 ' 43 ' 44 ' 45 Alter- 
natively, one can use a pure Gutzwiller formulation with 
an appropriate projector Pi^ The generalized energy 
functional for the Hubbard model reads: 

F(p, K ,D) = ^^(^jA^A,-*^) 

id 

i i 

which depends on the normal and anomalous parts of 
the density matrix and the variational double occupancy 
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parameters Di. In Eq. ([3]), pa = J2a Pif anc ^ we have 
introduced Nambu notation 



*l = (4 T , Ca ) *i= ( i 



-a 



(4) 



It is also convenient to define the following pseudospin 
vector 

Ji = 2*1 T ^ = 2 + Ci ^) ' (5) 

J? = ^*iT„*i - -| (4 T 4 4 - c a c lT ) , (6) 

J' = ^*!^*^^( c l T Cl T + c a c a-i), (7) 

where r$ denote the Pauli matrices. The components of 
Ji obey the standard commutation relations of a spin 

I 



1/2 algebra. We use boldface latter to indicate two- 
component vectors and 2x2 matrices in Nambu space, 
whereas we denote Cartesian vectors by using an arrow. 
The raising and lowering operators are defined as 

.it = jf+ijy = c\^ (8) 

jr = J; iJf r, /v . (9) 
Within these definitions the matrix A.j in Eq. ([3]) reads 



as 



A, = 



— 5 1 2 — JJ7) WT) w _ Z UJ 

2p^l^»T Z »U — 2 2 — ji~y 



and all expectation values (...) refer to the state 
The Gutzwiller renormalization factors are given by 



(10) 



y/Di - (Jf) - (Jj) y/1/2 + (J?) - A + *(S?) + y/1/2 + (Jf ) - Dj - a(Sf)y/Di - (J?) + (J. t ) 

v /l/4-((J,) + (^)) 2 



(11) 



r 



Note that in the limit (J^) = 0, where the matrix 
Aj is diagonal, one recovers the standard Gutzwiller 
energy functional as derived by Gebhard^ or Kotliar- 
Ruckenstein>iS 

F{p, D) = J2 UjaZiaZjaPji ~ P ^ P« + U Di > ^ 

i,j i i 

with the hopping renormalization factors 



_ V(Pu - A)(l - pa + Dj) + g - D l )D l 
VPu^-Pff) 

(13) 

The minimization of the energy functional of Eq. ([3]) 
leads to the stationary Gutzwiller wave function and to 
the associated stationary uncorrelated state \<fio). 



B. Calculation of pair fluctuations around general 
GA saddle points 

The energy functional of Eq. is a convenient start- 
ing point for the calculation of pair excitations on top 
of unrestricted Gutzwiller wave functions. Following the 
general approach of Refs. [l4lfl5lll6ll38l . we study the re- 
sponse of the system to an external time-dependent per- 
turbation which induces small-amplitude oscillations in 
the particle-particle channel: 



(14) 



Correspondingly, we have to expand the energy func- 
tional of Eq. (PJ around the stationary solution up to 



second order in the density and double-occupancy devia- 
tions. As already mentioned, we shall restrict to saddle- 
point solutions in the normal state: 



(Ji 



0. 



(15) 



and we remind the reader that a subscript or superscript 
indicates quantities evaluated in the stationary state 
\(/>o). Fluctuations are defined as 5p(t) = p(t) — po, etc. 

In the present case, particle-hole (ph) and particle- 
particle (pp) sectors in the expansion are decoupled and 
one obtains 



F = F + tr{h°5p} + SF ph + SF PP , 



(16) 



where we have introduced the Gutzwiller Hamilto- 
nian42^ 



Kf'lP,D] = ^S»,c 



(17) 



This coincides with the Kotliar-Ruckenstcin Hamiltonian 
matrix^ In particular, the diagonal elements (in the ba- 
sis of atomic orbitals) coincide with the Lagrange mul- 
tipliers of the Kotliar-Ruckenstein method, E^, after 
adding the chemical potential: 



Ei 



dF 
dpff 



P- 



(18) 



In Ref. l38| we have interpreted Ei CT as a local GA self- 
energy. 

Since we have included p in Eq. (TTJ) the eigenvalues of 
Eq. (|T7j> . £, describe the single-particle excitations with 
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respect to the chemical potential at the GA level. We de- 
note the particle (hole) energies above (below) the Fermi 
energy by £ pa (£_ ha ) with £ pa > > 

The transformation from real space fermions Cj CT to GA 
operators is written as 



Cia = £ 4>ia{jp)a po 



h 



4>ia(h)ahcr 



(19) 



and the amplitudes 4>ic{p) and <j)i a (K) correspond to the 
eigenf unctions of Eq. (|17|) and the index p (h) runs over 
empty (occupied) states. 

SF ph contains the expansion with respect to the 
double-occupancy parameters and the part of the density 
matrix, which commutes with the total particle number. 
This part of the RPA problem has already been stud- 
ied in detail in Refs. Ii~4lfl5l . where it was shown that the 
5D fluctuations can be eliminated by assuming that they 
adjust instantaneously to the evolution of the density ma- 
trix (antiadiabaticity condition). 

Finally, the particle-particle part of the expansion 
reads: 



SF pp = 



v^ki K^Kki 

ijkl 



(20) 



with Vijki = (Vkuj)* and the matrix elements of the ef- 
fective interaction are given by 



Van — 



5^ty {(<c] T c iT )o + (4t c jt)o) A° j++ A'- ++ 

3 

(<sW° + (cl lCjL )o) A<j__X!__} + ^ 



+ 1 ( c jl c il)0 

= - laA'iA'j |(c] T c lT ) + (c]i c a)o} for i^j 



Vijki 



\Vnij' 

= {Viiij 
= 



)* — t. tJ A a i++ A'j 



.A' 



for i =/= j 
for i j 
otherwise. 



of a superconducting ground state, one would have a cou- 
pling between (ph) and (pp) fluctuations and, therefore, 
the necessity to invoke the antiadiabaticity condition of 
Refs. [THUD to eliminate the SD deviations. 

Note also that in contrast to HF theory (where 
the expansion would be given by U^2 i 6{Ji)6{Jj~)), 
Eq. (|20|) contains correlations between distant pairs (i.e., 
5(J^)S(J~)) and also processes where pairs are created 

and annihilated on neighboring sites (i.e., <5(etj.ct^)). 

The remaining part of the formalism follows closely the 
particle-particle RPA as developed in nuclear physics (see 
Refs. l47ll48h . We define the v-th pp-RPA eigenstate of the 
N + 2 particle system by 



|JV- 
Rl 



■2,u 

p,p> 



Rt\N,0), 



VP 



Y hh' 

h,h' 



l h'\ 



,t 



(23) 
(24) 



where a pa and a ha create particles and holes in the single- 
particle levels of the Gutzwiller Hamiltonian Eq. (|17[) . 
The states \N,v) are unprojected states in the sense of 
Ref . flEI. i.e., they are auxiliary objects that have particle- 
particle RPA correlations but lack Gutzwiller correla- 
tions. This is because they result from creating particle- 
particle excitations on top of \cj>o) and not of |3>g). In 
the same way, the 77-th eigenstate of the N — 2 particle 
system can be represented as 



|Ar-2,J7)=i?t|iV-,0>, 



hh' 



i ii c 



t 

'h'\ 



pp' 



(25) 
(26) 



where \N, 0) is the unprojected RPA ground state of the 
-/V-particle system defined by 

R u \N,0) = R n \N,0)=0, 



Here, Ai TT i (with r, r' = ±) are the matrix elements of 
Eq. (jTUJ) and we have defined the following abbreviations 
for the derivatives: 



A' 



Ai T T — 



dA i+ _ 



dA, 



3^ Ai TT 



d(J+)d(J- 



(21) 
(22) 



where explicit expressions are given in Appendix [B] It 
is interesting to observe that, in contrast to the charge 
excitations in the particle-hole channel, the evaluation of 
the pair excitations can be performed without any as- 
sumption on the time evolution of SD. Only in the case 



and we adopt the convention that v (r\) runs over the 
n p (nfi) excitations of the N + 2 (N — 2) particle sys- 
tem. Within the pp-RPA scheme X and Y amplitudes 
can be associated with the following unprojected matrix 
elements: 



pp' 

I hh l 
A hh> 

y, 
pp 



(N,0\a p , iap1 \N + 2,v), 
-(N,0\a hl a h ^\N + 2,u) 
{N-2,r}\a h ^a hl \N,Q), 
-(N-2,r 1 \a p] a p , l \N,0). 



From the equations of motion for the amplitudes one can 
derive the following eigenvalue proble m 47 ' 48 
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Kpi 1 ^ ^P2t]^piP3^P2P4 VpiP2,P3P4 ^PlP2,/t3?i4 | ( Wp3P4 | _ [ \ ( Wj>lP2 j = 



(27) 



with J n the n x n identity matrix. The matrix ele- 
ments of the potential V can be derived from Eq. (j2T)|) 
by transforming to the GA representation with the help 
of Eq. dHJ). 

Equations (f2"T)) yield n,, + n/i eigenvectors which can 
be normalized as \W\ 2 — \Z\ 2 = ±1. The sign of the 
norm allows one to distinguish the n p addition eigenvec- 
tors (positive norm) from the n/j removal eigenvectors 
(negative norm) 4^ 

In the following, we will denote by E U (N) and E V (N) 
the A-particle energies of the Hamiltonian Eq. (fT]) when 
the fx term is absent and write the ll contribution ex- 
plicitly, so that the eigenvalues of Eq. (fT]) are given by 
E V (N) — llN . The eigenvalues and eigenvectors can be 
identified with the excitation energies and the amplitudes 
in the following form: 



E v (N + 2)-Eo{N)-2n, 



(28) 



for two-particle addition and 



E Q (N)-Er,(N-2)-2fi, 



(29) 



pp' pp r ' 

7 r > — Y n 



Stability requires that fi+ > > f2~ thus if there ex- 
ists a chemical potential ll such that these two conditions 
are satisfied the system is stable, otherwise a pairing in- 
stability arises. 

In a finite system, the allowed range of ll shrinks to zero 
as an instability is approached. From Eqs. (f2"8")) and (j2"5|) 
one sees that stability requires Eq(N + 2) — Eq(N) > 
E (N) - E (N - 2). This coincides with the stability 
against a transfer of a pair of particles from one clus- 
ter to another one in an ensemble. Notice that it does 
not necessarily coincide with the stability against single 
particle transfers. In the thermodynamic limit both con- 
ditions coincide with the well known stability condition 
that the compressibility must be positive. 

We are now in the position to evaluate the pairing cor- 
relation function within the TDGA. We are interested in 
on-site s-wave (s) pairing and intersite pairing: 



A 



1 



V2 



(30) 
(31) 



for two-particle removal. 



where S = x,y and i + S is a shorthand for the first 
neighbor of site i in the S direction. 

The dynamical pairing correlations can be computed 
from: 



3fM 



, , dte^<TA?(i)[A?(0)]n 



(32) 



1 „ (A,0|Af |A + 2. ,n [X + 2. i>\ i A / )< |..Y. 0} 



lu - m + io+ 



N. 



1 (JV,0|(Af)t|JV-2,7 ? )(JV-2,7 ? |Af|JV,0) 



U) — fir, 



i0+ 



where N s denotes the number of sites of the system and 
a = s,x,y. The average in the first line is done on the 
Gutzwiller projected RPA state. The latter is our best 
estimate for the true ground state of the system. The 
matrix elements in the second line are done on the un- 
projected states. In the spirit of the GA, Gutzwiller pro- 
jections are effectively taken into account by a renormal- 
ization of the operators. As usual local operators do not 
get renormalized, i.e., Af = Af, whereas non-local oper- 



ators acquire a renormalization through the z factors: 

Af = —f=(z i+ siz^ci + sic^ + Ziizi+s^ciia+s-]). (33) 

This is similar to the renormalization of the GA current 
operator in the computation of the optical conductivity J^. 
The validity of this prescription can be checked using sum 
rules as discussed below. Obviously in the HF theory 
based on BLA, bare operators are used in Eq. (f3"2")l . 

The matrix elements of the pairing operators can be 
obtained from the amplitudes X and Y using the trans- 
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formation Eq. ([19]) . In Eq. (|32|) the first and second term 
represent the correlations in case of two-particle addition 
and removal, respectively. 

In the case of a separable potential Vk 1 k 2 .k 3 k 4 . — 
Vk 1 k 2 v k 3 k 4 the pp-RPA equations can be easily solved. 
In particular for the case of a paramagnetic solution an- 
alyzed in Ref. H the TDGA interaction kernel Eq. (gOD 
only involves local pairs 



5F" = V'£ i 5{4,4 l )S{c i 



l c iV 



with 



V 



U-2H 
1-n 



(34) 



(35) 



where n = N/N s and S CT is defined in Eq. (fTSjl . In 
this case the on-site response can be found in momen- 
tum space: 

p £» = ^E^Wq^) ( 36 ) 

s q 

P aa (q,oj) = - — / die- t (TA«(t)[A«(Q)]t). 



with (a = s) 



q AT. ^ 



and is given by an effective ladder type equation: 

if.(q,«0 



f ss (q,w) 



i-ypp s ( q , w ) 



(37) 



(38) 



with /(£k) the Fermi distribution function, and r/k,k' = 
+ sign(£ k + £k')- The single particle energies are given 
by £k = £k — M where = z^e-^ + T,q is the GA dis- 
persion relation and e k = Y^j iij-e~ lk '( Rj ~ Ri *' is the bare 
one and zq is the hopping renormalization factor Eq. (113[) 
evaluated at the saddle point. 

We are interested also in the fluctuation response for 
intersite bond singlets with s-wave symmetry (a = +) or 
d-wave (a = — ) symmetry: 



Thus 

P aa (q,a;) - P QQ (q,^) + ^po^) ( 41 ) 

and the poles of P aa (q, u>) are determined by the local 
pair correlation (1 — FP° s (q, w) = 0) and the bare corre- 
lation. 

Eqs. ([3"5|) - (|4"Tj) are also valid in the BLA approach, re- 
placing the Gutzwiller local self-energy with the HF one 
(E ff = Un/2, z = 1) and taking V = U. 

The exact on-site s-wave spectral density satisfies the 
following sum rules: 



1 



dulm Pff (cj) = I - n + {n^rni) (42) 



dulmPs»(co) 



HI 



(43) 



In our case these sum rules allow us to compute RPA 
corrections to the double occupancy and they will be used 
below to evaluate RPA corrections to the GA ground 
state energy. Notice that the total spectral weight is 
equal to 1 — n. 

In addition the following energy weighted sum rule is 
satisfied 



Mf(oo) = - (T) G A+(2n-U)(l-n) (44) 



(45) 



with T the kinetic energy. This sum rule is satisfied ex- 
actly within the present pp-RPA in the sense that the 
right and the left hand side are equal, provided the ki- 
netic expectation value on the r.h.s. is computed at the 
GA level. Thus in contrast with Eqs. |g2J) and Q this 
sum rule provides no new information but is useful to per- 
form a consistency check. The same prescription is valid 
in the conventional HF plus pp-RPA approach, where 
the r.h.s. expectation values have to be computed at HF 
level4^ As shown in Appendix [C] the consistent renor- 
malization of intersite operators can be checked from an 
analogous sum rule. 



III. RESULTS 



A? = -(Af + Ar + aA» + aA-y) 

In this case the ladder series takes the following matrix 
form: 

£(q, w) = E° (q, u) + E° (q, «)E E(q, lj) (40) 

with 



P(n ,A- f P ssM Ps*(<l^)\ (vo 



In this section we present results for pair correlations 
in the repulsive Hubbard model within the framework de- 
veloped above. Since the RPA-type scheme used in the 
TDGA differs in some regard from the approaches usu- 
ally invoked in solid-state theory, we first illustrate the 
method for a two-site model which also can be solved ex- 
actly. Due to the mean-field character of the present ap- 
proximations the method is expected to improve with di- 
mensionality so this zero-dimensional example represents 
the worst case and allows us to give a first check of the 
potentialities and limitations of the method. Finally, we 
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compare our method with exact results on small clusters 
and demonstrate its superior performance as compared 
to the BLA. 



A. Two-site model 

In order to illustrate the formalism we consider a two- 
site model with two particles having up and down spins. 
The ground-state wave function reads as 



where 



|2,0) =a\S) + p\D) 



\S) = _(|1 T 2;)-|1;2 T )) 

\D) = -±=(|1 T 1;> + |2 T 2;» 



(46) 

(47) 
(48) 



and the corresponding amplitudes and the ground-state 
energy Eq(2) are given by 



2t 



E {2) 



-up 

U 
2 



1-1/1 



I6t 2 



(49) 
(50) 

(51) 



For four and zero particles there is only one state with 
energy E(A) — 2U and E(0) = 0, respectively. The en- 
ergy differences between two-particle addition (removal) 
states and the ground state are 



U 



0+ =£(2 + 2)-£ (2)-2 M = 
Q- = E (2) - E(2- 2) - 2[i = -Q+, 



16i 2 , 



and the chemical potential is taken as fj, = U/2 which is 
the exact value at half filling at an infinitesimal temper- 
ature. The corresponding matrix elements for local (s) 
and intersite (x) pairing operators 



A, 
A.,. 



1 

i 



(c u Ci T + c 2 xc 2T ), 
(c u c 2 t - Ci T c 2 x), 



(52) 
(53) 



read as 

<2,0|A s |2 + 2) = (2,0|At|2-2) = 0, (54) 
(2,0^2 + 2) = (2,0jA+|2 - 2) = a. (55) 

Notice that there is only one state with 2 ± 2 particles so 
we dropped the excitation index. One can check that the 
sum rules of Eq. (j42|) -(|44l are satisfied. For example the 
double occupancy is given by (3 2 /2 and the first moment 
sum rule 



2 

— IT) 



(2fj,-U)(l-n) 



n+\(2, 0|A S |2- 
2{U-u a )[3 2 = 



2)| 2 - 
-<T) 



fi-|(2,0|At|2-2>| 5 



is also fulfilled, as it should. 

Consider now the same model within the TDGA. On 
the GA level, one finds two single particle states for each 
spin direction which at half filing can be put as: 



iha = ~t(l - U 2 ) + S CT - (J,, 



i(l-u 2 ) + £ c 



(56) 
(57) 



where the h-states is occupied with a spin-up and spin- 
down particle and we have defined u = U /Ubr with 
Ubr = 8t. For U < Ubr there is a paramagnetic so- 
lution which becomes insulating at the Brinkmann-Ricc 
transition point Ubr^ For u > \[2— 1 the more stable 
solution is an antiferromagnetic broken symmetry solu- 
tion which does not have a Brinkmann-Rice transition 
point. 

For the paramagnetic solution one has spin- 
independent Lagrange multipliers which from the 
Kotliar-Ruckenstain (or Gebhard's) scheme are obtained 
as 



v - v - U 



(58) 



Since there is only one two-particle addition and removal 
state, we have to diagonalize the following 2x2 pair 
fluctuation RPA problem 



2i(l - u 2 ) + V/2 V/2 

V/2 2t(l - u 2 ) + V/2 



= n 



1 

-1 



(59) 



r 



and the interaction V corresponds to the local part in the expansion Eq. (|20|) (note that the first derivatives of 



8 



A® vanish for a paramagnetic solution) 
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V = -2t(l - u 2 )A" + 



=Uu(2-u)l^ (60) 
1 — n 1 — u 



where A" is defined in Appendix [Bl Here, the diverging 
part, proportional to 1/(1 — n), is canceled by an analo- 
gous contribution in the first term. The same result can 
be obtained from Eq. ([35]) by taking the limit n — >• l, 38 
and coincides in modulus with the effective interaction 
at half filling in the particle-hole case. 13 ' 14 ! 15 This is con- 
sistent with the fact that the attractive-repulsive trans- 
formation converts particle-particle fluctuations in parti- 
cle hole-fluctuations with equal interaction but sign re- 
versed^ 

Diagonalization of Eq. ([59|) yields the two eigenvalues 



fi ± =±2t(l-u 2 ). 



V 



2t(l-u 2 )' 



For the eigenvectors, the following relations hold 

2t(l-u 2 ) 



{X ± +Y ± ) 2 = ±- 



±\2 



{X±Y _ (y±) 



±i, 



(61) 



(62) 
(63) 



from which one can compute the amplitudes for the lo- 
cal and intersite pairing correlations between unprojected 
states: 

(2,0|A S |2±2) = -L(X±+y±), (64) 

(2,0|A JB |2±2> = -±(X ± -Y ± ). (65) 

Finally, the expectation values between Gutzwiller pro- 
jected states are the same as for unprojected states in 
the case of A s and should be renormalized in the inter- 
site case, i.e., Aa, = (1 — u 2 )A x . See Appendix ICl for a 
consistency check of this prescription. 

For the local pairing operator we find that the first 
moment sum rule Eq. (|44D 

n+|(2,0|A s |2 + 2)| 2 -fT|(2,0|At|2-2)| 2 (66) 

= 2i(l - u 2 ) = ~{T)ga - 2/i(l - n) + U(l - n) 

is satisfied (note that n — 1) in the TDGA, as antici- 
pated. 

The TDGA two-particle addition and removal ener- 
gies are displayed in Fig. [T] and compared with the exact 
ones. Solutions obtained in the paramagnetic regime are 
shown with solid lines. Remarkably, due to a cancella- 
tion of (1 — u) terms the Brinkmann-Rice transition does 
not reflect in the excitation energies of the paramagnetic 
phase, which become soft at a much larger value of the re- 
pulsion (u — 1 + \/2) • Instead, the transition shows up in 
the TDGA pairing correlations (right panel of Fig.[l|. In- 
deed, the Gutzwiller renormalization factors Zi a drive the 
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FIG. 1: (color online) Left panel: Excitation energies for two- 
particle addition and removal computed for the two-site model 
with two particles. Right panel: Matrix element of the local 
(lower curves) and intersite (upper curves) pairing operator. 
The underlying GA saddle-points are paramagnetic (p) (solid 
lines) and antiferromagnetic (dashed lines). The critical value 
u — v2— 1 of the corresponding transition is indicated by the 
vertical dotted line. 



matrix elements of the intersite pairing operator to zero 
at the Brinkmann-Rice transition point. It is interesting 
that neglecting zi a in the pairing operator, the suppres- 
sion is replaced by a divergence so that the Zi a cancel the 
unphysical divergence but overcorrect it. This problem 
is partially cured in the broken symmetry state (dashed 
lines) but for U — > oo the matrix element tends to 1/2 
whereas the exact result approaches | (2 — 2| A.^ |2, 0} | 2 — 5> 1 
since the ground state can be written as Aj.|2 — 2) in this 
limit. Instead, the ground state with broken symmetry 
has only one configuration of the two that generate the 
singlet ground state which explains the 1/2 factor. For 
the local pairing operator the behavior is better. 

We notice that the TDGA excitations energies are in 
quite good agreement with the exact values, especially 
when one allows for the broken symmetry solutions. 

Finally we can use Eq. (|42|) to compute the pair fluc- 
tuation derived double occupancy: 

D RPA = (n iT n a ) = (2, 0|c uCiT |2 + 2)(2 + 214^^2,0). 

(67) 

Fig. [5] (left panel) shows the double occupancy in dif- 
ferent approximations. HF completely neglects correla- 
tions and hence the double occupancy is independent of 
U . The BLA based correction drives the approximation 
much closer to the exact result at small U but strongly 
underestimates the correction at large U. In the GA, 
correlations are taken into account already at the static 
level and the double occupancy is strongly suppressed as 
a function of U . Thus the pp-TDGA correction is small 
and brings the double occupancy close to the exact one 
in a much larger range of interaction. The fact that the 
pp-TDGA is close to the GA double occupancy indicates 
that the theory is nearly self-consistent. In Ref. l38jit was 
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shown that this feature is enhanced in two-dimensions 
pointing to an improvement of the performance as the di- 
mensionality is increased. In contrast the BLA is clearly 
quite far from being self-consistent in this sense. 

For large U, the exact double occupancy is of order 
UJ ~ t 2 /U 2 due to the same charge fluctuations that 
build the double exchange interaction J Since the GA 
and the TDGA results in Fig. [2] are for paramagnetic so- 
lutions the double occupancy vanishes at the Brinkman- 
Rice point. This makes clear that the scale J is not 
present in the paramagnetic GA or TDGA. 

The above results allow us to compute the pp-RPA 
ground state energy using the coupling constant integra- 
tion trick^^ 



Eq PA = —2t+ UN S I dU'D RPA (U'). (68) 



Here the first term is the ground state energy for U = 0. 
Restricting on the paramagnetic solutions, we find for the 
TDGA and BLA 



TDGA 



BLA 



= -2t + At{yJl + 2u-u 2 -1), (69) 
= -2t + 2t(y/l + U/2-l), (70) 



which both yield the exact result luq up to second order 
in U/t. 
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FIG. 2: (color online) Double occupancy (left panel) and 
ground-state energy (right panel) of the two-site Hubbard 
model computed within HF, BLA, TDGA, and exact diago- 
nalization. Only paramagnetic ground states are considered. 

Due to the more accurate estimate for the double oc- 
cupancy (see Fig. fj^) it turns out that the TDGA gives 
a significantly better approximation for the ground-state 
energy (see Fig. [2b) than the BLA over the whole range 
of U, before the onset of the Brinkmann-Rice transition. 
This should be compared with an analogous calculation 
in the particle-hole channel in Ref. [l5|. In the latter case 
one obtains a singularity in E at u — %{y2— 1) due to 
the onset of the antiferromagnetic ground state. Since 
there is no instability in the particle-particle channel, 



such problems do not arise in the present case. Thus it 
seems convenient in general to compute RPA corrections 
to the energy in a channel free from instabilities. 



B. Comparison with exact diagonalization 

In this section, we study the dynamical pairing corre- 
lations within the TDGA on small clusters and compare 
our results with the BLA approach and exact diagonal- 
izations. 

We start by computing the long-distance pairing corre- 
lations (Af (A|)t) for bond singlets. Within the pp-RPA 
these correlations are obtained integrating the addition 
part of the dynamical spectral function: 



(Af (Af 



1 



Thus the comparison with the exact diagonalization re- 
sults provides a stringent test of the total spectral weight. 
We remind the reader that in the TDGA the singlet op- 
erator is renormalized by the factors. 

In Fig. [5] we show results for 10 particles on a cluster 
with 18 sites, tilted by 45°, and, therefore, having all the 
spatial symmetries of the infinite lattice. A particular 
representation of this cluster, with its boundary condi- 
tions, is shown in the upper panel of Fig. [3J In this case, 
the GA ground state is paramagnetic. We concentrate on 
non-overlapping singlets and the distance between them 
is measured from their centers. 

The two singlet operators can form a perpendicular 
configuration, like 's-b' in the upper panel of Fig. 02 or a 
parallel configuration, like 's-a' in the same figure. Notice 
that for Rij = Ri — Rj = 2 there are two points in the 
lower panel corresponding to two parallel configurations 
in which one of the singlets is displaced either along the 
x- or along the y-direction (labeled 'a' in the upper panel 
of Fig. p. 

In Fig. O the vertical bars point to the GA value 
of the singlet correlations, i.e., the decoupled result of 
(A| (Ajy)aA but still renormalized with the Zi a factors. 
In the same spirit of Ref. |22j the length and orientation 
of the bars reflects the 'vertex contribution'. This quan- 
tity measures the correlation induced interaction between 
two singlets which is attractive when the TDGA value is 
larger than the one computed within the bare GA. For 
nearby singlets we observe excellent agreement between 
the exact diagonalization result and the TDGA for both 
U/t = 4 and 10 (see lower panels of Fig. [3]). In this case 
we also observe an effective attractive interaction. In gen- 
eral, with increasing distance TDGA overestimates the 
exact correlations, and the difference becomes more pro- 
nounced for larger U/t. This behavior can be expected 
due to the fact that the Gutzwiller method is based on a 
local projector which neglects intersite correlations. One 
can therefore anticipate that the incorporation of inter- 
site projections (as in Jastrow-type wave function o 55 i 56 ) 
into the TDGA would lead to an improvement of the 
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long-distance pair correlations. Nevertheless, it turns out 
that the TDGA yields a rather good description of long- 
distance pair correlations (especially for moderate values 
of U/t) as compared to exact diagonalizations. 
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FIG. 3: (color online) Top panel: Sketch of the 18-site clus- 
ter and corresponding boundary conditions. Bottom panels: 
Bond singlet pairing correlations on the 18-site lattice for 
U/t = 4 and 10. The bars at the TDGA symbols point to 
the bare GA value of the singlet correlations. The separa- 
tion between bond singlets is defined as the shortest distance 
between their centers. Taking the singlet 's' as reference, it 
should be noted that there exist two parallel singlet configu- 
rations 'a' with distance R s — R a — 2. The singlets 'b' and 
's' form a perpendicular configuration. 

In order to compare dynamical properties of the TDGA 
with BLA and exact results, we compute the d-wave and 
s-wave correlations of Sec. Ill Bl The comparison of the re- 
sults in different approximations can be done by aligning 
the chemical potentials (as in the last figure of Ref. 1381) or 
aligning the absolute energies. In the following, we adopt 
the last procedure by eliminating the chemical potential 
from the response function. We define u>' = uj+2h so that 
the poles in Pg a (u') occur at J = E V {N + 2) - E (N) 
and J = E V {N - 2) - E (N). 

In Fig. 21 we show the addition spectra for P^ a (u') 
evaluated for a 18-site cluster with 10 particles and 
U/t = 10. This corresponds to a closed shell configura- 
tion where the HF and GA approach yield paramagnetic 
solutions. The on-site s-wave case P^ s {lo') has already 
been analyzed^ for a smaller cluster and it is shown in 
Fig. 2] for 18 sites a as reference since, as explained in 
Sec. IIIB1 it also determines the overall features of the 
intersite pairing correlations. Fig. [5] (main panel, dashed 
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FIG. 4: (color online) Imaginary part of the (two-particle ad- 
dition) pairing correlation function of Eq. (|32p for on-site pair- 
ing (a), extended s-wave (b) and d-wave (c) symmetry. Re- 
sults are for 10 particles on a 18-sites cluster and U/t = 10 for 
which the underlying mean-field solution in BLA and TDGA 
is paramagnetic. The delta-like excitations are convoluted by 
lorentzians with width e = 0.2t. The arrows in the upper 
panel point to the chemical potentials obtained within the 
various methods. The insets detail the frequency evolution of 
the weight at small energies. 



dotted line) displays the effective interaction V for the 
present case as a function of U / 1. 

In Ref. [38| we have shown that the on-site pair excita- 
tions for large U in both BLA and TDGA are dominated 
by an antibound which lower band edge is at ui' = U . 
The band states at energies It ~ 5t have very small 
spectral weight. In the top panel of Fig. 0] the full 
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black arrow indicates the value of 2/i in the exact case 
with fi = (fi+ + n~)/2 and fi+ = E{N + 1) - E(N), 
fj,- = E(N) - E(N - 1) . The position of 2/i in the GA 
(red dashed arrow) is very close to the exact one whereas 
in HF (blue dotted arrow) it is shifted to higher energies. 
Thus, aligning the chemical potentials, the position of the 
BLA (TDGA) antibound state is in disagreement (agree- 
ment) with the exact result^ In addition, the TDGA 
gives a good account of the low-energy spectral weight 
(see inset), which is much larger and shifted to higher 
energies in the BLA. This was understood as due to the 
different way the self-energy is renormalized by interac- 
tions in GA and in HF42> This dramatic difference in per- 
formance gets greatly amplified for intersite correlations 
[see Fig. \Mjo) and (c)] because the band states acquire 
significant weight. 




u/t 



FIG. 5: (color online) Main panel: Local part of the effective 
interaction Vu : a (see Eq. ([20])) for the 18-site cluster. The 
solid and dashed lines refer to the half-filled cluster with para- 
magnetic and spin-density wave ground state, respectively. 
The dashed dotted line is for the 10-particle system with a 
paramagnetic ground state. Upper-left inset: Effective inter- 
action between pairs on nearest-neighbor sites Vajj. Lower- 
right inset: Effective interaction Vu,ij between a local pair 
on site Ri and an intersite pair on the bond defined by Ri 
and Rj. Vu.ij > (<)0 for Sf < (>)0 (solid and dashed lines 
respectively) . 

The intersite correlations have also significant weight 
at the energy of the s-wave antibound state w' ~ U. This 
can be understood from the second term of Eq. (HT|) and 
corresponds to processes in which an intersite pair decays 
into an on-site pair. Notice that the d-wave pair cannot 
mix with a local s-wave pair at q = but couples at finite 

q 

In addition, the exact result shows a satellite at oj' « 
2U for the intersite cases corresponding to a process in 
which the two particles are created at neighboring sites, 
that are initially single occupied, thus leading to two on- 
site pairs in the final state. This can be visualized as the 
creation of the two particles in the upper Hubbard band. 
Since this band is not present in our starting point (i.e., 
GA), this satellite is absent in the pp-RPA. This feature 
reappears if one starts from a GA state that has both 



lower and upper Hubbard bands at the cost of sponta- 
neous symmetry breaking. This case is analyzed in the 
following. 

We now consider the half-filled system where both HF 
and GA have a spin-density wave (SDW) ground state. 
The dynamical pairing correlations, for 18 particles in 
the 18-site cluster, are shown in Fig. [6] Again, we are 
comparing absolute energies, but the position of 2/i in the 
different approximations is aligned because the chemical 
potentials coincide and are equal to the exact result /i = 
U/2. 

For local s-wave, extended s-save, and d-wave fluctua- 
tions the exact result has a broad distribution of weight 
centered around u>' ~ 2U. Clearly this corresponds to 
the high-energy satellite of the previous case. Now prac- 
tically all sites are singly occupied so the probability to 
find an empty site where to create a local s-wave pair 
is very small and there is practically no weight at the 
energy of the antibound state ui' ~ U. Indeed, if the 
TDGA computation is done with a paramagnetic state 
(which cannot reproduce the satellite at 2U) the response 
becomes identically zero/22, This can be also seen from 
the sum rule of Eqs. (|42[) . since in the paramag- 
netic TDGA and GA the double occupancy vanishes for 
the half-filled system above the Brinkman-Rlce transition 
point. 

Breaking the symmetry and allowing for the SDW, 
both BLA and TDGA give a good estimate of the energy 
scale but with a more narrow distribution of weight. No- 
tice that the on-site s-wave response has a much smaller 
intensity consistent with the fact that the feature origi- 
nates from intersite excitations which get mixed with the 
on-site response. 

Processes in which the two particles are created on dif- 
ferent singly occupied sites are allowed in both approx- 
imate theories and they lead to approximately similar 
results. In this case the interaction is practically ineffec- 
tive. In fact, the excitation spectra evaluated from P° 
(from the bare Gutzwiller Hamiltonian) shown with thin 
lines in Fig. [5] differ only slightly from the TDGA result 
and the latter only shows a small shift of spectral weight 
to higher energies due to the inclusion of particle-particle 
correlations. This is in contrast to the paramagnetic case 
away form half filling (see Fig. [J) where all high-energy 
(antibound) states are determined from poles in P(u>') 
but are absent in the non- interacting case P (u) 1 ) as dis- 
cussed before. 

The insets of Fig. [5] show the evolution of the two- 
particle spectral weight. For both extended s-wave and 
d-wave symmetry the total weight approaches 1/2 in the 
pp-RPA theories in the limit U — > oo as for the AF so- 
lution of the two-site model in the right panel of Fig. [1] 
This can be easily understood from the Neel limit. In 
contrast, the exact result has a larger weight which can 
be understood from quantum fluctuations that induce 
spin flips. In the extreme case of the two site model, the 
fluctuations lead to a local singlet and the exact spectral 
weight for U — > oo is twice the approximate one as dis- 
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FIG. 6: (color online) Imaginary part of the (two-particle ad- 
dition) pairing correlation function Eq. (|32p for on-site pairing 
(a), extended s-wave (b) and d-wave (c) symmetry. Results 
are for a half-filled 18-site cluster and U/t — 10. The under- 
lying mean-field solution in BLA and TDGA is a spin-density 
wave. The delta-like excitations are convoluted by lorentzians 
with width e = 0.5i. The insets detail the frequency evolution 
of the two-particle spectral weight. 



cussed in connection with Fig. [T] Clearly this result is 
closely related to the reduction of the magnetization in a 
quantum antifcrromagnet. 

The structure of the TDGA interaction kernel Eq. p7|) 
is more complex than the paramagnetic case. Fig. [5] 
displays the corresponding non-vanishing elements as a 
function oiU/t. For small on-site repulsion the local con- 
tribution behaves as Vu t a ~ U , independent of the filling 
and the ground state, as it should. With the onset of the 



SDW at U/t w 4.1, Vu t ii (dashed line) starts to deviate 
from the corresponding interaction in the paramagnetic 
limit (full line) which diverges at the Brinkmann-Ricc 
transition. 38 Interestingly, the interaction between local 
pairs on adjacent sites (upper left inset of Fig. [5J is al- 
ways attractive in the SDW phase with a maximum at- 
traction at U/t m 10. In addition one finds an attractive 
or repulsive interaction between a local pair on site i and 
an intersite pair Ci^Cj.^. For i and j nearest neighbors 
the interaction is attractive (repulsive) if the pair has the 
same (opposite) spin with respect to the underlying Neel 
magnetic moments of sites i and j. Nevertheless, these 
additional fluctuations in the TDGA cannot overcome 
the strong residual on-site repulsion so that the system 
remains stable against a transition towards superconduc- 
tivity. This also holds for a SDW ground state away from 
half filling. 

Finally, as in the two-site example, we calculate the 
energy correction from the TDGA for a half filled 4x4 
cluster from Eqs. (|4*2"|) and In Fig. [71 we compare the 
corresponding result for the particle-particle and particle- 
hole channel (from Ref. ll4T ) with the bare GA and the ex- 
act ground-state energy. The underlying saddle-point of 
the GA solution is a SDW. As can be seen from Fig.[7J the 
particle-hole TDGA correction is approximately twice 
that in the particle-particle channel. The former gives a 
quite accurate approximation for intermediate values of 
the on-site repulsion but tends to overshot the exact re- 
sult for large U/t (note that these energy corrections are 
not derived from a variational principle and thus do not 
constitute an upper bound for the exact result). With the 
considered range of U/t the particle-particle corrections 
always are slightly higher in energy than the exact ground 
state. However, in comparison with the HF+RPA energy 
corrections which are by far too large?^ the TDGA yields 
a reasonable approximation to Eq in both particle-hole 
and particle-particle channel. 



IV. CONCLUSIONS 

In this paper, we have extended the TDGA towards 
the inclusion of pair correlations in the Hubbard model. 
The present analysis is complementary to previous com- 
putations in the particle-hole channe l 14 ! 15 ' 16 where we 
have analyzed the spectrum of charge- and spin excita- 
tions. In comparison with exact diagonalization results 
the TDGA yields a very good agreement for the dynami- 
cal pair correlation function, especially for the low-energy 
excitations and away from half filling where it performs 
significantly better than the HF based BLA theory. 

Compared to numerical methods 5 -- our approach can 
be pushed to much larger systems. In particular, it is 
suitable for the evaluation of pair correlations in the 
negative-?/ Hubbard model, where the approach is at 
least qualitatively capable to capture the crossover from 
weak coupling BCS to strong coupling behavior^ An 
important outcome of this investigation concerns the jus- 
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FIG. 7: (color online) Ground state energy of the half-filled 
4x4 Hubbard model computed within GA, TDGA in the 
particle-hole channel, TDGA in the particle-particle channel, 
and exact diagonalization. 



tification of the antiadiabatic assumption for the time- 
evolution of the double occupancy, which was needed in 
the charge particle-hole channel. 

For the Hubbard Hamiltonian with particle-hole sym- 
metry (e.g., nearest- neighbor hopping at half filling), the 
superconducting instability and the charge-density wave 
instability in the particle-hole channel are degenerate. 
Evaluation of the latter instability within the TDGA re- 
quires to invoke the antiadiabaticity condition whereas 
the expansion in the particle-particle channel of Eq. (j20|) 
goes without it. The fact that the two calculations within 
the TDGA for the particle-hole and the particle-particle 
sector give the correct degeneracy of the instabilities for a 
charge-rotational invariant system clearly indicates that 
the antiadiabaticity assumption was indeed correct, i.e., 
other possibilities, as keeping the double occupancy fixed 
at the stationary value (rather than to follow the time 
evolution of the density matrix), would have led to an 
unphysical breaking of charge-rotational symmetry. 

Our approach does not produce a superconducting in- 
stability in the Hubbard model. This is because, in the 
paramagnetic phase, the effective interaction does not 
contain the attractive part due to the exchange of spin 
fluctuations. Indeed, the superexchange scale J is absent 
in the paramagnetic GA. On the other hand, it turned 
out that the effective interaction in the SDW phase con- 
tained attractive contributions between nearest-neighbor 
pairs which, however, are too weak in order to yield a 
superconducting instability once the SDW is present. 

The flexibility of the present approach allows one to 
study collective modes in inhomogeneous superconduct- 
ing states once the Hamiltonian is augmented with a suit- 



able pairing potential. This may find application in the 
physics of high-T c cuprates, where in many compounds 
the occurrence of electronic inhomogeneities (e.g., in the 
form of stripes) is now well established. Another possi- 
ble application is in the field of ultra-cold atoms where 
the ground state is intrinsically inhomogeneous due to 
the presence of the confining parabolic potential. In ad- 
dition, recent studies incorporate disorder to produce a 
glassy state, which may lead to similar physics as in the 
cuprates 
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APPENDIX A: DERIVATION OF THE 
CHARGE-ROTATIONALLY INVARIANT GA 
WITHIN THE SLAVE-BOSON APPROACH 

The procedure implemented in the following essentially 
consists of three steps. Assume that in our initial refer- 
ence frame we have non-vanishing superconducting or- 
der, i.e., (J,j + ) 7^ and (or) in (Jr) ^ 0. First, we rotate 
them locally to a new frame where these expectation val- 
ues vanish, i.e., (Jj ) = {J7) = 0. This allows, as a 
second step, the introduction of slave-bosons and associ- 
ated fermions within the Kotliar-Ruckenstein scheme. 
For the bosons, we apply the saddle-point (mean-field) 
approximation. Finally, in a third step, we rotate the 
fermions back to the original reference frame. 

We define the local rotations in charge space by the 
following transformations 

= U}* 4 *J = 9f\Ui, (Al) 

where 

Uj = cos(^i/2)l + ism(ipi/2)T ■ 17, (A2) 

and 77 = (r) x ,T) y ,0) is the rotation axis of length unity. 
The inverse transformation reads as 

$i=Ui«i *t = *tut. (A3) 

Within the first step of our procedure we have the 
requirement that the transformed pseudospin vector is 

given by Jj = (0, 0, J*). Applying the transformation of 
Eq. (|A3|) to this vector one obtains the following relations 

Jf = -%sin(^)Jf, (A4) 

Jf = Vx M<Pi)Ji*, (A5) 

Jf = cos(^)Jf. (A6) 
Note that the spin operator 

S? = \ - l) = S? (A7) 
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is unchanged by the transformation. 

Since by definition off-diagonal order vanishes in the 
rotated frame we can now, as a second step, apply the 
Kotliar-Ruckcnstcin slave-boson scheme: 



C-itr %i<7 fia 



t 



with 



(A8) 
(A9) 



4 e i + Pl-aPi-o 



The double (d), singly (p CT ), and empty (e^) occupancy 
bosons are constrained by the following relations: 



25? = 2S? 



(A10) 

(All) 
(A12) 



Since we follow essentially a Gutzwiller-type approach, 
we now apply the mean-field approximation for the 
boso ns. With the help of Eqs. flSjj), ([A10]i . (|ATT|) 
and (|A12|) we can eliminate all bosons but di and ex- 
press them via expectation values in the original refer- 
ence frame and c&. One finds 



(J? 



PI 



cos(^) 

1 1 

2 cos(ipi) 



[J? 



(A13) 
(A14) 



and a = ±1 in the latter equation. 

Summarizing the steps we have performed so far the 
original fermion operators Ci a are related to the Kotliar- 
Ruckenstein transformed ones fi a in the rotated frame 
via the transformation 



*i = 



J 



= W 




(A15) 



with 
W = 



z iiT cos^ i(r] x — ir] y ) sin ^j.j. 



(A16) 

Finally we transform the fermion operators fi a back to 
the original frame [see Eq. (|A1|) ] 




U 



/it 



(A17) 



so that the charge-rotational invariant Gutzwiller repre- 
sentation of the fermions is given by 



fi r 
f 



wu 




(A18) 



The complete transformation matrix A = WU reads as: 



with 




tanV-^ 



(J? 



(A20) 



and J+, J i and J? are defined as in Eqs. ([7]), fl5|), and © 
but with the pseudofermion operators fi a instead of Ci a . 

Note the formal similarity of Eq. (JTUJ) with the cor- 
responding transformation of the spin-rotation invariant 
Gutzwiller approximation^ 

We now turn to the transformation of the interaction 
term of the Hubbard Hamiltonian which has to be per- 
formed within the second step of the above scheme but 
before the saddle-point approximation has been applied 
to the bosons. Rewriting the interaction in terms of the 
transformation Eq. (|A15|) and obeying the constraints 
Eqs. (|A10 )) -((A12 |) we obtain 



n i\ n %i = cos 2 -£-d\di + sin 2 -^-e]ei, 



2 fi t 



(A21) 



which now can be again treated in the mean-field approx- 
imation. 

We finally obtain for the charge-rotational invariant 
Gutzwiller energy functional of the Hubbard model 



+ U Yj d 2 i - J t z (^Jl + ta,n 2 <pi-l) . (A22) 



The term multiplying U is clearly the sum of the double 
occupancies Di of the original hamiltonian. Therefore, 
we perform the substitution 



Di 



tan 2 (fi 



(A23) 



to reach Eq. ([3]). With this definition it is also straight- 
forward to proof the equivalence of the saddle-point z- 
factors Eq. (|A9[) with the Gutzwiller renormalization fac- 
tors Zi„ from Eq. (fTTj) and of the two representations of 
the transformation matrix A in Eq. (fT0|) and Eq. (|A19|) . 



APPENDIX B: DERIVATIVES OF THE 
HOPPING FACTOR 



The derivatives appearing in Eqs. (|2"Tj) and 

jiven by: 



are 
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A' 



dAi+- 
d{J~) 




%TT 



A - d(J+)d(J, 



o 
1 



2Jf 



(1 - rii) 2 Hi — 1 n iT (l - m T ) 



1 - rn y/n iT (l - n iT ) 



Dj 



i - m + A 



1 /'I - II, -<r D, 



TLj T Di 



(Bl) 
(B2) 
(B3) 



Notice that the index r for the matrix element of A appears at the place of the spin index on the right. In this case 
one should interpret + =] and — ={. 

In case of a homogeneous, paramagnetic saddle point these expressions simplify to 



A' = 



A- 



2z 



n(2 — n) 1 — n 



l-n + D+ ^JD{\ -n + D) n-2D 



(B4) 
(B5) 



where zq denotes the spatially homogeneous Gutzwiller 
factor defined in Eq. (flB"]) . 



APPENDIX C: SUM RULE FOR INTERSITE 
PAIRING OPERATORS 

We define local and extended s-wave pairing operators 
on a hypercubic lattice (coordination number Z) as 

A r = [^i c jt + c Ji c d 

3 

where 7^ = 1 for ij nearest neighbors and 7^ = oth- 
erwise. We consider the Hubbard model Eq. ([I]) with 
nearest-neighbor hopping only. The kinetic energy oper- 
ator can then be represented as 



the following sum rule 



(CI) 



and the commutator of A' with H yields 

[H,Al\ = tAf -UA\. (C2) 
For the double-commutator with (A^ s )'!' we obtain 
[(Af)t, [H, A',]] = -2Zt + Zt{n A + n lL ) 

Et t 

nm 



(C3) 



which only depends on one-particle operators. 

Taking the expectation value of Eq. (|C3[) and inserting 
a complete set of N + 2 and N —2 particle states we find 



£<[(Ar)t,[/f,A!]]) 
^-(^\(Ar)^- 2 )(^- 2 \^\K) 



= -2N L Zt{l-(n)) + j(T) 

(C4) 

where Nl denotes the number of lattice sites, (n) the 
particle density and lj± = ±[E V {N ± 2) - E (N)]. 

In case of the time-dependent GA the states |\l/^ r±2 ) 
are interpreted as Gutzwiller projected RPA states, i.e. 
they incorporate Gutzwiller and particle-particle correla- 
tions. In the spirit of the GA matrix elements are evalu- 
ated in term of the unprojected RPA states as 



«|(Af)t|^ 



N-2\ 



(0,7V|(Af)V,iV-2) 



where A^ s is renormalized with the z factors as in 
Eq. (|33|) . As in HF+RPA theories one expect that the 
r.h.s is equal to the l.h.s when evaluated at the static 
mean field level in this case the GA. Indeed evaluating 
the kinetic expectation values on the GA approximation 
(incorporating the z factors) one finds and identity as 
expected. 

We can explicitely demonstrate the procedure for the 
half-filled two-site model where local and intersite pairing 
operators are defined in Eqs. (|S"2"|) and (|S"3")) . For this 
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special case the sum rule Eq. (|C4[) becomes Evaluating the r.h.s. of Eq. (IC5[) within the GA one 

obtains 

{[(A es )\[H,A 1 ]]) = w+<*g|A , |* a+2 )(* 2+2 |(A M ) t |*g) 

- "-(•8KA-)t|»»><^|AK) ^^(-^(l-^-^l-^) (07) 

= j t (T). (05) 

For the l.h.s. of Eq. E3 wc find from Eqs. El, and El so that in the time-dependent GA the intersite pairing 

operator Eq. (|53[) has to be renormalized according to 

u + {0, 2\A l \2 + 2) (2 + 2|(A") t |0, 2) 

A K = -Wl - u 2 )(c u c 2T - ci T c 2i ) (C8) 



c-(0,2|(A-)t|2-2>(2-2|A i |0,2) x _ 1 M . 



= - V [(a+) 2 - (y+) 2 ] + i w - [(x-) 2 - {Y-f} ^ 

= -- [w+ + w~] = -2S = -J7. (06) in order to fulfill the sum rule Eq. ([CSjl . 



1 M.C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963). 

2 G. Baskaran, Z. Zou, and P. W. Anderson, Solid State Com- 
mun. 63, 973 (1987). 

3 G. Kotliar and J. Liu, Phys. Rev. B 38, 5142 (1988). 

4 A. Paramekanti, M. Randeria, and N. Trivedi, Phys. Rev. 
Lett. 87, 217002 (2001). 

5 S. Sorella, G.B. Martins, F. Becca, C. Gazza, L. Capriotti, 
A. Parola, and E. Dagotto, Phys. Rev. Lett. 88, 117002 
(2002). 

6 W. Metzner and D. Vollhardt, Phys. Rev. Lett. 59, 121 

(1987) . 

7 W. Metzner and D. Vollhardt, Phys. Rev. B 37, 7382 

(1988) . 

8 F. Gebhard, Phys. Rev. B 41, 9452 (1990). 

9 M.C. Gutzwiller, Phys. Rev. 134, A 923 (1964); ibid. 137, 
A 1726 (1965). 

10 G. Kotliar and A.E. Ruckenstein, Phys. Rev. Lett. 57, 1362 
(1986). 

11 J. Lorenzana and G. Seibold, Phys. Rev. Lett. 89, 136401 

(2002) . 

12 J. Bunemann, F. Gebhard, T. Ohm, S. Weiser, and W. 
Weber, in Frontiers in Magnetic Materials, ed. A. Narlikar, 
Springer (2005). 

13 D. Vollhardt, Rev. Mod. Phys. 56, 99 (1984). 

14 G. Seibold and J. Lorenzana, Phys. Rev. Lett. 86, 2605 
(2001). 

15 G. Seibold, F. Becca, and J. Lorenzana, Phys. Rev. B 67, 
085108 (2003). 

16 G. Seibold, P. Rubin, F. Becca, and J. Lorenzana, Phys. 
Rev. B 69, 155113 (2004). 

17 J. Lorenzana and G. Seibold, Phys. Rev. Lett. 90, 066404 

(2003) . 

18 G. Seibold and J. Lorenzana, Phys. Rev. Lett. 94, 107006 

(2005) . 

19 G. Seibold and J. Lorenzana, Phys. Rev. B 73, 144515 

(2006) . 

20 E. Dagotto, J. Riera, and A. P. Young, Phys. Rev. B 42, 
2347 (1990). 

21 F. Becca, A. Parola, and S. Sorella, Phys. Rev. B 61, 
16287(R) (2000). 

22 S.R. White, D.J. Scalapino, R.L. Sugar, N.E. Bickers, and 
R.T. Scalettar, Phys. Rev. B 39, R839 (1989). 

23 A. Moreo, Phys. Rev. B 45, 5059 (1992). 



24 N. Bulut, D.J. Scalapino, and S.R. White, Phys. Rev. B 
47, 14599 (1993). 

25 S. Zhang, J. Carlson, and J.E. Gubernatis, Phys. Rev. 
Lett. 78, 4486 (1997). 

26 E. Plekhanov, F. Becca, and S. Sorella, Phys. Rev. B 71, 
064511 (2005). 

27 R. Arita and K. Held, Phys. Rev. B 73, 064515 (2006). 

28 T. Yanagisawa, Phys. Rev. B 75, 224503 (2007). 

29 RA. Ferrell, J. Low Temp. Phys. 1, 423 (1969). 

30 D.J. Scalapino, Phys. Rev. Lett. 24, 1052 (1970). 

31 H. Takayama, Prog. Theor. Phys. 46, 1 (1971). 

32 S.R. Shenoy and PA. Lee, Phys. Rev. B 10, 2744 (1974). 

33 M. Cini, Solid State Commun. 20, 605 (1976). 

34 GA. Sawatzky, Phys. Rev. Lett. 39, 504 (1977). 

35 M. Cini and C. Verdozzi, Solid State Commun. 57, 657 
(1986). 

36 For a review, see C. Verdozzi, M. Cini, and A. Marini, J. 
Electron Spectrosc. Rel. Phenom. 117, 41 (2001). 

37 K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. 
Hecker Denschlag, A.J. Daley, A. Kantian, H.P. Biichler, 
and P. Zoller, Nature 441, 853 (2006). 

38 G. Seibold, F. Becca, and J. Lorenzana, Phys. Rev. Lett. 
100, 016405 (2008). 

39 V. M. Galitskii, Zh. Eksp. Teor. Fiz. 34, 14 (1958); [Sov. 
Phys. JETP 7, 104 (1958)]. 

40 J. Kanamori, Prog. Theor. Phys. 30, 275 (1963). 

41 R. Fresard and P. Wolfle, Int. J. Mod. Phys. B 6, 237 
(1992). 

42 B.R. Bulka, Phys. Stat. Sol. B 180, 401 (1993). 

43 B.R. Bulka and S. Robaszkiewcz, Phys. Rev. B 54, 13138 
(1996) 

44 M. Bak and R. Micnas, J. Phys.: Condens. Matter 10, 
9029 (1998). 

45 J.O. Sofo and GA. Balseiro, Phys. Rev. 45, 377 (1992). 

46 J. Bunemann, F. Gebhard, K. Radnoci, and P. Fazekas, J. 
Phys. Condens. Matter 17, 3807 (2005). 

47 P. Ring and P. Schuck, The nuclear many-body problem 
Springer- Verlag, New York, 1980. 

48 J. Blaizot, G. Ripka Quantum theory of finite syatems, 
MIT Press, 1986. 

49 D.J. Thouless, Nucl. Phys. 22, 78 (1961). 

50 W.F. Brinkman and T.M. Rice, Phys. Rev. B 2, 4302 
(1970). 



17 



R. Micnas, J. Ranninger, and S. Robaszkiewicz, Rev. Mod. 
Phys. 62, 113 (1990). 

J. Lorenzana, G. Seibold, and R. Coldea, Phys. Rev. B 72, 
224511 (2005). 

G.D. Mahan, Many Particle Physics, Plenum, New York, 
1990. 

A. L. Fetter and J. D. Walecka, Quantum theory of Many 

Particle Systems, McGraw-Hill, New York, 1971. 

R. Jastrow, Phys. Rev. 98, 1479 (1955). 

Recently, the role of the long-range Jastrow factor for hav- 



ing an accurate variational wave function in the Hubbard 
model has been discussed by M. Capello, F. Becca, M. 
Fabrizio, S. Sorella, and E. Tosatti, Phys. Rev. Lett. 94, 
026406 (2005). 

E. Dagotto, Rev. Mod. Phys. 66, 763 (1994). 

F. Giinther and G. Seibold, Physica C 460-462, 1041 
(2007). 

L. Fallani, J.E. Lye, V. Guarrera, C. Fort, and M. Inguscio, 
Phys. Rev. Lett. 98, 130404 (2007). 



